function [X, L, X_err, L_err, iter, flag] = ...
    aitkenpower(A, Xo, X_ref, L_ref, maxit, tol, print)
%
%  [X,L,X_err,L_err,iter,flag] = power(A, Xo, X_ref, L_ref, maxit, tol)
%
%  SORPOWER applies the Aitken-accelerated power method to the eigenvalue 
%  problem:
%     
%    (A - L*I)*X = 0
%
%  The unaccelerated (power method) iterates take the form
%    X[0]   = normalized guess, i.e. ||X[0]||_2 = 1.0
%    Z[l+1] = A * X[l]                                   (1)
%    L[l+1] = Z[l+1]' * Z[l+1]                           (2)
%    X[l+1] = Z[l+1] / sqrt(L[l+1])                      (3)
%
%  Aitken's delta-squared process gives a better estimate for
%  x2 in a sequence x0, x1, and x2.  Here, 
%
%  Acceleration replaces (3) by
%    X[l+1] = X[l] + w[l] * ( Z[l+1] / sqrt(L[l+1]) - X[l] )
%
%  A w[l] of 1 yields no acceleration.  A w[l] >(<) leads to
%  over(under)relaxation.  Note, 0 <= w <= 2.
%
%  Reference:
%
%    E. Z. Mueller
%    Laban-PEL documentation
%
%  Inputs:
%
%    real       A(N,N)     -- N-by-N matrix
%    real       xo(N)      -- Initial guess N-vector.
%    real       lambda_ref -- The actual eigenvalue for checking errors
%    integer    maxit      -- The maximum number of iterations.
%    real       tol        -- The error tolerance (defined...?)
%
%  Outputs:
%
%    real       x(N)       -- The computed eigenvector.
%    real       lambda     -- The computed eigenvalue.
%    real       rho        -- The estimated dominance ratio.
%    real       err_norm   -- The norm of the error.
%    integer    iter       -- The number of iterations performed.
%    integer    flag       -- flag (0=okay, 1=maxit reached, -1=bad bad news)
%

if (nargin==6) 
  print = 1;
end

flag    = 1;
X_err   = 0;
L_err   = 0;

% Normalize the initial guess.
Xo      = Xo / norm(Xo, 2);

% Generate first three.
[X0, L0] = powermethod(A, Xo, X_ref, L_ref, 1, tol, 0);
[X1, L1] = powermethod(A, X0, X_ref, L_ref, 1, tol, 0);
[X2, L2] = powermethod(A, X1, X_ref, L_ref, 1, tol, 0);

% Set initial extrapolation parameters.
w1 = 1;
w2 = aitken(L0,L1,L2);

for iter = 3:maxit
    % Save.
    L0 = L1;  
    L1 = L2;
    w1 = w2;
    % Update.
    Z    = A * X2;         
    L2   = norm(Z, 2);
    w2   = aitkenp(L0,L1,L2);
    if ( check_w(w2, w1) == 1 )
        if (print==1) 
            disp(sprintf(' ** Using w = %6.5e at iter %5i ',w2, iter));
        end
        % We use w2 to update L2 and X2
        L2 = L1 + w2 * (L2 - L1);
        X2 = (1 - w2) * X2 + (w2 / L2) * Z;
    else
        % No overrelaxation
        X2   = Z / L2;
    end
    L = L2;
    X = X2;
    X_err = norm(X - sign(X(1)/X_ref(1))*X_ref, 2);
    L_err = abs(L - L_ref);
    if ( X_err <= tol && L_err <= tol)
        flag = 0;
        break;
    end;
    if ( mod(iter, 300) == 0 )
        printout(iter,X_err,L,L_err,print,abs(1-w2/w1));
    end
end

% final messages and results
if (flag==0 && print == 1)
    disp(' *** PI: Final Results *** ');
end
if (flag==1 && print==1)
    disp(' *** Warning: PI did not converge ***');
end
printout(iter,X_err,L,L_err,print,abs(1-w2/w1))

end

function printout(iter,X_err,L,L_err,print,wrat)
if (print==1)
  disp(sprintf(...
  ' Iter = %5i, X_err = %6.5e, L = %6.4e, L_err = %6.3e, abs(1-w2/w1) = %6.4e ', ...
  iter, X_err,L,L_err,wrat))
end
end

% checks for "exponential mode"
function y = check_w(w2, w1)
    if ( abs(1-w2/w1) < .5 )
        y = 1;
    else
        y = 0;
    end
end

% vectorized aitken extrapolant
function x = aitkenp(x0, x1, x2)
    x = (x1-x0)./(2*x1 - x2 - x0);
end
